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We introduce a new method to reconstruct unknown quantum states out of incomplete and noisy 
information. The method is a linear convex optimization problem, therefore with a unique minimum, 
which can be efficiently solved with Semidefinite Programs. Numerical simulations indicate that 
the estimated state does not overestimate purity, and neither the expectation value of optimal 
entanglement witnesses. The convergence properties of the method are similar to compressed sensing 
approaches, in the sense that, in order to reconstruct low rank states, it needs just a fraction of the 
effort correspondig to an informationally complete measurement. 
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C\) . I. INTRODUCTION 

All we can know about a quantum system can be compactly represented by its state or density matrix, i.e., a trace 
O one positive operator (/?). The observable properties of the system are represented by Hcrmitian operators (O), and 
their expectation values are given by the rule: Tr(pO). Consider the system is defined on a Hilbert space of dimension 
d. In this case, the state (p) is a d x d Hermitian matrix, or a real d 2 dimensional vector in the Hilbert-Schmidt 
space. As the state is normalized (Tr(p) = 1), it has just d 2 — 1 independent real parameters, so it is in fact ad 2 -l 
I"** , dimensional real vector. Therefore, a carefully chosen set of d 2 — 1 independent measurements, in the sense that it 
i i ' forms a complete basis in the Hilbert-Schmidt space, allows for the reconstruction of the state. This approach is 
usually referred to as quantum tomography. 

There are many methods to estimate a quantum state (for a review, see, for example, [1]). An important class of 
these methods is known as maximum-likelihood quantum tomography (MaxLik) (sec, for example, [2-5]). MaxLik is a 
statistical inference approach which yields probability distributions as close as possible to the measured frequencies. 
Nonetheless, these probability distributions frequently correspond to non-positive operators, due to the always present 
experimental imprecisions. It is important to stress that our method does not follow this line. Our heuristic consists 
of performing a variational search in the state space for the best trace one positive operator compatible with the 
measured data. Note also that MaxLik cannot guarantee that fake entanglement be not attributed to the system 
[6]. Our numerical simulations indicate that the method we propose does not overestimate entanglement. Another 
7-H . important aspect of our method is that it can be solved exactly, for it can be cast as a linear convex optimization 
problem, while the usual algorithms for quantum state tomography, though convex, are non-linear, and the actual 
algorithms implementing them can be plagued with many local minima. 

To perform the complete set of measurements necessary to reconstruct the quantum state of a moderately sized 
system could be impossible in practice. As an example, consider the experimental generation of the 8-qubit entangled 
state reported by Haffner et al. [9]. To tomograph the 8-qubit system, 65536 projectors needed to be measured (2 16 ), 
and the computational effort to apply MaxLik to such a state is discouraging, though it was done. But the tomography 
of a larger system could be inaccessible experimentally, and yet we would like to assess its properties. The exponential 
growth of the Hilbert space is insurmountable, in this respect. Using the Maximum Entropy Inference Principle of 
Jaynes [7] (MaxEnt), one can estimate a quantum state out of incomplete information. In this approach, one searches 
for the maximum entropy state, compatible with the measured data [8] . This state corresponds to maximum ignorance 
with respect to the unmeasured set of observables forming the complete basis for the density operator. The Horodecki 
[6] showed an example where the direct application of MaxEnt yielded fake entanglement, i.e., the reconstructed state 
was entangled, while the available data were compatible with a non-entangled state. Our method can also estimate a 
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state out of incomplete information, and with the advantage of not overstimating entanglement, as indicated by our 
simulations. 

Another interesting aspect of our method is that it can identify incompatible data. Suppose a set of complete 
measurements are being performed, aiming a quantum state reconstruction. It can happen that, for a particular 
complete measurement, some parameters of the experimental setup have fluctuated to a point of changing the quantum 
state, but it went unnoticed by the experimentalist. In this case, data for different states have been collected. When 
our method is applied to these measurements, the algorithm will detect that the data are incompatible with a single 
quantum state. Then some of the data can be discarded, and the state can be reconstructed. 

In the next section, we develop a variational method to estimate a state, which relies on incomplete and noisy 
information of the quantum state. The method we obtain has the particular form of a linear convex optimization 
problem, known as Semidchnitc Program (SDP), for which efficient and stable algorithms are available [10-12]. 
Though this unfortunate name can suggest a kind of black-box computational algorithm, a SDP is not so. It consists 
of minimizing a linear objective under a linear matrix inequality constraint, precisely, 

minimize c'x 

{m 
F(x) = F + X i F i ^ °> W 
i=l 

where c £ C m and the Hcrmitian matrices Fi £ C nxn are given, and x £ C m is the vector of optimization variables. 
F{x) > means that F(x) is a positive matrix. The problem defined in Eq.l has no local minima. When the 
unique minimum of this problem cannot be found analytically, one can resort to powerful algorithms that return 
the exact answer [11]. To solve the problem in Eq.l could be compared to finding the eigenvalues of a Hcrmitian 
matrix. If the matrix is small enough or has very high symmetry, one can easily determine its eigenvalues, but in 
other cases some numerical algorithm is needed. Anyway, one never doubts that the eigenvalues of such a matrix can 
be determined exactly. We point out that we have successfully used SDPs before, in the development of powerful 
methods to construct entanglement witnesses [13-16]. 

After discussing the methodology, we present a section with numerical examples that we consider representative, 
namely, a full rank highly symmetric two-qutrit density matrix, and a pure (rank one) five-qubit state. Our method 
works equally well in these two extreme cases. Then we show reconstruction of four-qubit mixed states of all ranks, 
in order to have a better taste of the convergence properties of our method. These numerical examples suggest that 
our method converges very fast for low rank states, or for states that, though of high rank, have very high symmetry. 
These convergence properties are similar to recent introduced approaches in compressed sensing [17]. To reconstruct, 
with high fidelity, full rank states with no symmetry, an informationally complete measurement is needed. But even in 
these difficult cases, our method can yield reasonable lower bounds to entanglement out of incomplete information [16]. 
We illustrate this feature by evaluating optimal entanglement witnesses [13-16] in qubit-qutrit full rank states, with 
increasing number of measurements. Besides the numerical examples, we have also successfully tested our approach 
in a real experiment, where entangled qutrits were generated in a quantum optics setup [18]. 



II. THEORY 



The state of a quantum system is represented by its density matrix p, which is a d 2 — 1 dimensional real vector in 
the Hilbert-Schmidt space. Therefore, we can write it as: 

d+ld-l d 2 -l 

p=c i+j2J2 CiiP i =c ° I+ J2 c ^' ( 2 ) 

1 = 1 i=l A=l 

where the P\ = Pa arc Hermitian operators, forming a complete basis in the Hilbert-Schmidt space, / is the identity, 
and A is just a convenient compact index. 

Now, for concreteness, we assume the operators in each class I are rank one projectors, forming a complete mea- 
surement, i.e., 

(3) 
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d 2 — 1 

Co, which is not an independent parameter, in this case is given by Co = (1 — Sa=i Ca)/cL Note that, in Eq.2, we 
use just d—1 projectors from each complete measurement labelled by I (for the d probabilities sum to one). Besides 
being convenient for our calculations, this basis also has the minimum number of projectors necessary to expand a 
Hermitian matrix with fixed trace, but the method we derive is not basis dependent. For instance, if one chooses to 
expand the state in a SU basis, like a generalized Bloch representation, 

d 2 -l 

P=- L I+ H n<Tl ' W 

i=l 

to obtain the coefficients j-j = Tr(pa i ) 1 one should consider the spectral decomposition cr, = Ylj=i a y ana ^ we 

are back to the measurement of projectors, Tr(\ij)(ij\p). But now, instead of the d 2 — 1 projectors of Eq.2, one has 
d x (d 2 — 1) projectors, i.e., one complete measurement for each operator <n. The use of POVMs (Positive Operator 
Valued Measure) also poses no difficulties, for it consists of measurements of rank one projectors again. 

We assume that from a set of d 2 — 1 projectors (or observables), forming a complete basis in the Hilbert-Schmidt 
space (the space of linear operators), just N(< d 2 — 1) have been measured. We refer to measured and unmeasured 
projectors (observables) as the known and unknown sets, respectively. Now we introduce the following cost operator 
(the sum of the projectors in the unknown set), which can be thought of as a Hamiltonian: 

d 2 -i 

H= ^ P A . (5) 

We want a trace one positive operator, which minimizes the cost function: 

d 2 -i 

E = Tr(Hg) = ]T q X , (6) 

where q\ are positive numbers corresponding to the projections of g on the unknown set. Our variational principle 
now reads: 

SE = Tr(HSg) = 0. (7) 

The Lagrangian for this problem reads: 

N 



L= E + a[Tr(g)-l] + J2MTr(gPx)-px] + 

(8) 

\=N+1 

where a, (3\, and 7a are Lagrange multipliers corresponding to the constraints on the quantum state. 
Eq.2 in matrix form reads 

P = C T P, (9) 

where C is a column vector with the real coefficients C\, and P a column vector collecting the matrices Pa- Consider 
the overlap matrix S, with elements S^ v = Tr{P^P u ), and the column vector V, with elements p\ — Tr{pP\). Then 
we have 

C = S'- 1 P. (10) 

Only when the probabilities in V are exact, the vector C yields a positive operator, by means of Eq.9. But it is never 
the case, in practice. The frequencies obtained from an experiment are noisy, thus: 

Tr(pP x ) E[px-e, Px +e}, (11) 



with e positive, and hopefully small. 
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To account for Eq.ll in our algorithm, we introduce the additional constraints: 

Tr(gP x ) > (1 - A\)p\, 

Tr(gP x ) < (l + A x ) Px , (12) 
A A > 0, A G [l,N]. 

Again, we want to determine a trace one positive operator g, with minimum A\. It is important to note that Eq.12 is 
not to be interpreted as any kind of statistical treatment. By minimizing A\, we are simply adjusting the frequencies 
obtained in the laboratory, such that we get a positive matrix by means of Eqs.9 and 10. This works only when Eq.ll 
is valid for all frequencies. If some frequencies came from another state, say Tr{p other P x ) 1 then the algorithm will 
diverge, and we have identified incompatible data. To build an appropriate Lagrangian for this problem, we introduce 
positive numbers v x and w\, which should be minimized, and rewrite Eq.12 as: 

Tr(gP x ) - (1 - A\)p\ = v x , 

(1 + A\)p\ - Tr(gP x ) = w x . ^ 

The Lagrangian now reads: 

N 

L = E + J2(&x+v x +w x ) + a[Tr(g) - 1] + 

A=l 
d 2 -l 

lx[Tr(gP x ) - q x } + 



X=N+1 
N 



(14) 



^{Ca[(1 + A\)p\ - Tr(gP x ) - w x } + 



Vx[Tr(gP x )-(l-A x )p x -v x ]}, 

where ( x and r] X are additional Lagrange multipliers. Note that this Lagrangian defines a convex problem. We want to 
minimize a linear function, constrained by matrix inequalities. Our variable is q, and the search space is restricted to 
the cone of positive matrices. This is a typical convex optimization problem, known as Semidcfinite Program (SDP), 
as described in the Introduction (Eq.l). SDPs have a unique minimum, and can be solved very efficiently [10]. Our 
SDP reads: 

minimize (E + J2 X =i ^a) 
Tr(g) = 1, 

subject to{ A x > 0, (15) 
Tr(gP x ) > (1 — A x )p x , 
Tr(gPx) < (l + A x ) Px , VA e [1,7V]. 

Now we have concluded the derivation of our method. Eq.l 5 returns a state g which is the optimal approximation to 
the unknown state p. 



III. APPLICATIONS 



Now we would like to illustrate the use of Eq.15. As the examples show, our method does not overestimate purity 
and neither the expectation value of optimal entanglement witnesses [13-16], and it can also identify incompatible 
data. The examples also show that few measurements are needed to reconstruct both low rank states, and full rank 
states with high symmetry (Figs. 1, 2 and 3). In Figs. 1, 2 and 3, the set of measurements we chose for the 
state expansion is mutually unbiased bases (MUB) [20, 21], in the sense that any two vectors of different complete 
measurements have the same overlap's absolute value. As demonstrated by Wootters et. al [21] and Ivanovic [20], this 
is the best informationally complete projective measurement one can do. It is optimal both in the statistical sense and 
in the number of projectors to be measured. The POVMs that would be equivalent to MUBs arc the Symmetrically 
Informationally Complete POVMs (SIC-POVM) [22]. But while the MUBs are known for all Hilbcrt spaces which 
have dimension of power of a prime, SIC-POVMs, though conjectured to exist in all dimensions, are known just in a 
few particular cases. In Fig. 4, we used the SU(6) basis, to illustrate that our method works with any kind of basis. 
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FIG. 1: Estimation of some properties and reconstruction of a two-qutrit full rank Werner state (/3 = —0.8, cf. Eq.16). We plot 

Purity, Tr(g 2 ), Fidelity, Tr(\J pi gpi), Trace Distance, ^Tr\g — p\, and entanglement, Tr(WgQ~) (Wg is the trace-one optimal 
entanglement witness), against the number of measured projectors. The measurements have statistical errors of up to 50%, 
according to a uniform distribution. In the last two panels, we introduced measurements from another Werner state (/3 = +0.8), 
to show that the algorithm is capable to identify incompatible data (see text for details). 



Our first example is a highly mixed 2-qutrit Werner state [19] (purity=0.23), which is full rank, and has entanglement 
-0.21, according to its optimal trace one entanglement witness [13, 14]. The explicit expression of this state is: 

p w (16) 

with — 1 < /3 < 1. pis separable for /3 > — |. F is a swap operator for two qutrits, 

F= E (17) 
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FIG. 2: Reconstruction of a five-qubit low entangled pure state. We plot Purity, TV (g 2 ), Fidelity, Tr(Vphp^), Trace Distance, 
\Tr\g — p\, and entanglement, Tr(WgQ) (Wg is the trace-one optimal entanglement witness), against the number of measured 
projectors. Note that the exact state is recovered with just 160 measurements out of the total of 1056, that would be necessary 
to completely span the Hilbert-Schmidt space. 



For two qutrits, we have ten complete measurements, with nine projectors in each one. In Fig.l, we plot, against 
the number of measured projectors, purity, fidelity and trace distance to the true or exact state, and the witnessed 
entanglement {Tr(WgQ) - the more negative is this expectation value, the more entangled is the state). In the 
calculations of entanglement, we use an optimal trace one witness, obtained with techniques developed in [13-16]. The 
optimal entanglement witness has to be calculated for each state individually. The (simulated) measured frequencies 
(Eq.ll) have statistical errors of up to 50%, according to a random uniform distribution. In the first panel, one 
can see that with about twenty measurements, from a total of ninety, we have a practically perfect estimation 
of the entanglement, and the state was reconstructed with high fidelity, or low trace distance. Note that all the 
properties we are calculating converge monotonically, never been overestimated. In the second panel, the third 
complete measurement (indexed 19 to 27) corresponds to a different Werner state. Note that we had convergence 
up to the 18 t/l projection, and then all the properties started to diverge, to resume convergence again in the fourth 
measurement. In the last panel, the incompatible data was moved to positions 81 to 90, and one can see convergence 
for the measurements before this. The fluctuations seen in all the panels are due to the large statistical errors we 
imposed, added to the fact that without all the complete measurements, there is always the possibility of a family of 
close states to fit the data. 

To highlight the efficiency of our method, in Fig. 2 we applied it to a a five-qubit low entangled pure state (it is a 
normalized random vector with 32 complex coefficients). In principle, one needs to perform 33 complete measurements 
(1056 projectors) to reconstruct such a state, but we needed just five (160 projectors). The entanglement plotted in 
the figure is the genuine five-party one, given by the trace one optimal entanglement witness [13-16]. The tomography 
calculation took six seconds in a laptop, with a 2.8GHz processor, and 3GB of RAM, running MATLAB under Linux. 
The five-party entanglement calculation was more expensive, and took about forty minutes. 

Fig. 3 illustrates the convergence properties of our method, as the rank of the density matrices increases. For a 
system of four qubits (Hilbert-Schmidt space dimension of 256), we consider one hundred random density matrices 
of each rank, varying from pure states (rank one) to full rank states (rank 16). We plot the average number of 
measurements to reconstruct the state with high fidelity (Tr\p — g\ < 10 ), against the rank. In the figure we also 
indicate the minimum (Best Case) and maximum (Worst Case) number of measurements needed in the reconstruction, 
in each sample. We see that the state reconstruction needs very few measurements for low rank states, and it can 
need all the measurements for high rank states. This behavior is similar to the compressed sensing approach reported 
in [17]. 

As a last example, we consider a qubit-qutrit system. The states will be expanded in the 36 SU(6) observables (<7j) 
forming an informationally complete basis in the Hilbert-Schmidt space (cf. Eq.4). The Hamiltonian (cf. Eq.5) to be 
used in Eq.15 is formed by the eigenprojectors of the 36— N unmeasured <7j (cf. Eq.4). For 10 3 random full-rank density 
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FIG. 3: Reconstruction of random four-qubit entangled states of all ranks. The average number of measurements needed to a 
highly faithful reconstruction (Tr|p— g\ < 10~ 6 ) is plotted against the rank of the state. For each rank, it is employed a sample 
of one hundred states. For each sample, it is also plotted the minimum (Best Case) and maximum (Worst Case) number of 
measurements needed. 
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FIG. 4: Reconstruction of 10 3 full rank random qubit-qutrit entangled states. We plot the fraction of the true entanglement, 
on average, as the number of measurements increases (see text for details). The average trace distance to the true state is also 
plotted. We see that an almost perfect reconstruction is possible only with all the measurements, but non null lower bounds 
to the entanglement are obtained with less measurements. 



matrices (pi), we plot, against the number of measured observables, the average trace distance (10~ 3 X)i=i Tr\g~i — 
Pi\/2), and average fraction of entanglement, defined as follows. The entanglement of a state is E(p) = Tr(W p p), 
if this expectation value is negative, and zero otherwise. W p is the optimal trace one entanglement witness for the 
state p, calculated with the techniques described in [13-16]. Thus, the Entanglement plotted in the figure is given by 
10 -3 5^«=i E(Qi) I E(pi). The full rank states will need all the 36 measurements to be reconstructed with high fidelity, 
for they have no symmetry. Even though, Fig. 4 shows that a non null lower bound to the entanglement is obtained 
with less measurements, and of course this bound tends to the correct entanglement as the number of measurements 
increases. 



IV. CONCLUSION 



We developed a method that yields an estimate (g) of a unknown quantum state (p), without the need to perform 
an informationally complete measurement. Our numerical simulations indicate that our method does not overesti- 
mate entanglemnt, and does not underestimate purity. These quantities tend to the true values, as the number of 
measurements increases towards an informationally complete measurement. Low rank states, or high rank states with 
high symmetry can be reconstructed with high fidelity, using few measurements. This is simply because the number 
of independent parameters in such density matrices is much less than d 2 , the dimension of the Hilbert-Schmidt space. 
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Note that it is a line of investigation we started in [16], in the context of entanglement detection with few measure- 
ments. The method can be useful in the study of larger systems, where an informationally complete measurement is 
out of question. In this respect, we note that the convergence properties of our method is similar to the compressed 
sensing approach recently introduced in [17]. 

We conclude by mentioning that our method has been successfully employed in a real experiment reported in [18]. 
We also mention that we have extended our method to the problem of process tomography [23]. 
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